SIGNAL
PROCESSING & SIMULATION NEWSLETTER
Discrete Fourier Transform (DFT) and the FFT
Let’s
take this continuos signal composed of three sinusoids.
Assume
that fa = 1, fb = 3, and fc = 5. The waveform
is plotted below.
Figure 1 - A signal representing a square wave
We
can draw the Fourier transform of this signal easily by examining the
amplitudes of each of these frequencies and then putting one-half on each side
of the y-axis as shown in Fig. 2 for a two sided spectrum. Real signals such as
this one produce only one sided spectrums also shown below.
Figure 2 - The two-sided and single-sided Fourier Transform of g(t)
This
is the theoretical Fourier Transform of the continuos waveform g(t). The
Fourier Transform tells us that there are just three frequencies in the signal
and no others. There is no ambiguity in the results. This ideal Fourier
Transform is what we want to see when do the Fourier Transform on a analyzer or
on a computer but in reality this is nearly impossible to obtain. All
implementations of the Fourier Transform are attempts to achieve the
theoretical results, however, digital signal processing introduces
approximations and truncation effects which keep us from realizing the ideal.
Figure 3 shows the outline of the same signal along
with dots that represents what we actually see of the signal on a oscilloscope.
This is because most signals we capture are sampled versions of the real analog
signal. We pulse the analog signal every so often and then plot these sampled
values. We connect the samples and get a proxy to the signal.
Figure 3 - The discrete samples of a real signal as shown by dots. We
really do not know the underlying shape of the signal.
Mathematically, the sampled signal is obtained by
multiplying the target signal with an impulse train of period t. Since we usually collect only a limited
number of samples we limit the length of the impulse train to a certain time
window. The discrete signal is expressed as
So before we can even look at a signal, two things
have happened. 1. we have multiplied the target signal by an impulse train and
made the continuos signal a discrete signal and 2. we have chosen to collect only
a limited number of samples, in effect windowing the sequence with a
rectangular window function.
Fig 4a shows the original signal and its Fourier
Transform. In (b) we have the Fourier Transform of a pulse train which is used
for sampling the original signal. The Fourier transform of the impulse train
consists of just one frequency, the sampling frequency.
Next step is the rectangular window that limits the
infinite impulse train. Its Fourier Transform is shown in c and is the
well-known sinc function.
Now we have a signal which is a product of three
signals.
for t < T
1 2 3
The first is the original signal, the second is the
impulse train of period t and the third is a step
function lasting for time T. What about the Fourier Transform of the product of
these three signals? Mathematically we know that multiplication in time domain
of two signals results in convolution of their spectrums in the frequency
domain. So we have an inkling that the convolution of all three of the Fourier
Transforms may not give us the spectrum in (a). But is it close enough to (a),
and if not how different is it?
Figure 4 - The quandary of digital signal processing
a. The continuos target signal, b. is multiplied by an impulse train of
frequency fs, c. The limiting window of length T, d. The sampled and
limited signal, what do we get here?
Let’s
define some terms.
Sampling frequency: fs
Sampling frequency is a measure of how often we
pulse the continuos signal to obtain the samples. The quantity sampled is the
amplitude of the signal.
Sample Time: t
Time between each sample. It is also the
inverse of the sampling frequency.
If sampling frequency, fs = 20
samples/sec, then
t
= 1/fs = .05 sec
Total Number of Samples: N
The total number of samples collected or
available = N
Sequence Time Length: T
This the time length of the
collected sequence and is equal to the product of the sample time and the total
number of samples.
= N t
Sample index: k
The
sequence time is no longer continuos so instead of t, we use a discrete time
measure called kth sample. This is an index of the samples. Its
range extends from 0 to N-1, where N is the last sample. Each kth
sample of total N samples, is located at time k times t secs.
Kth
Figure 5 - Referring to individual samples
In a continuos signal we refer to a particular point
at its instantaneous value of t. For discrete signals, We refer to any
particular sample as g(kt). So each sample differs in time from the previous one by t secs. For example the 3rd
sample similarly is located at (3 x .05) = .15 secs and the 10th
sample is located at time (10 x .05) = .5 secs..
Define the g(t) in its discrete form as
g(t) = g(k t)
Harmonic index: n
Frequencies that are integer
multiple of a fundamental frequency are referred to as harmonics of that
fundamental frequency. In computing DFT, we use the concepts of harmonics in a
special way.
From the sampling theorem, we know
if we want to recreate a signal from its samples then we must collect at least
twice its frequency number of samples per second. This also says that we can
detect frequencies in a signal only up to one half of its sampling frequency.
So the values of n ranges from
This says that we can only detect
half as many harmonics as the total number of samples. However the index itself
goes from -(N-1) to +(N-1) and spans both sides of the spectrum reflecting the
positive and negative components of the frequency.
Relating the Discrete
Fourier transform to the Continuous Fourier Transform
Recall
that the Fourier transform is given by
-1-
The
above equation says that if we multiply the target signal by a complex sinusoid
of harmonic frequencies 1, 2, 3, ..n
one at a time and then integrate the results, the integration yields the
amplitude of the nth harmonic. Why? Because multiplication by the
sinusoid acts as filter for all other frequencies. (Refer to Fourier Transform
Tutorials No. 1 and 2). The objective is to compute the complex coefficients cn,
which when plotted give us the frequency content of the signal.
= amplitude of the nth
harmonic sine wave
= amplitude of the nth
harmonic cosine wave
the complex
coefficients which produce the spectrum
The
integration limits in the Fourier transform formula of Eq. 1 go from - in to +
in. What does that
mean for our sampled sequence of N samples?
Fourier
Transform also requires that the signal be periodic. But looking at only N
samples we can not tell if the samples cover one exact period, more than one or
less than one period. In order to do the Fourier Transform, we need at least
one whole period or the result is suspect.
We
already see some problems as we go from continuous to discrete processing, The
problems are
1.
We
do not have an infinitely long series and
2.
We
do not know if the N samples we have observed cover a single period,
less
than a single period or more than one.
Let’s
continue despite the fact that we don’t know if what we are about to do is
right. We are going to assume that these two things will not cause us much
trouble and the results will be acceptable.
Now
let’s change time from continuos to discrete by making the following
substitutions for time and frequency.
g(t) = g(k t)
Equation 1 becomes
-2-
We
have also changed the integral to summation to note the change from continuous
to discrete as two processes are equivalent in the two domains.
Fundamental frequency of the
signal
By
applying the Fourier Transform algorithm on these N samples, we have made an
implicit assumption. We have assumed that the signal is periodic over N
samples, so we have assumed that the fundamental frequency of our signal is
equal to the inverse of time T of the N samples. We express the fundamental
frequency as
We
can rewrite T as a function of sample time, t and total number of samples
chosen, N to alternately express the fundamental frequency in terms of fs
and N.
This
is a very important concept to understand. It says that you have artificially
set the fundamental frequency to the sampling frequency divided by the total number
of samples observed. It is a strange idea seemingly having nothing to do with
the target signal and in fact this is true.
This
frequency referred to as the fundamental
frequency of the signal really is kind of a resolution frequency and has
nothing to do with the target signal. It just means that we resolve the target
signal components in integer multiples of this resolution frequency.
Let’s
say that we sampled the above signal at sample time of 1/20th sec
and observed 60 samples. Then the fundamental frequency is
Now
when we compute the Fourier Transform we will be stepping this fundamental
frequency by integer multiples. With f0 = .333, the next harmonic
would be f1 = .666 and so on. The harmonics used in the analysis are
not integers but integer multiple of
the fundamental frequency of the signal as determined by the sampling frequency
and the N samples observed. An alternate way to see these harmonics is see them
as bins which collect energy. In DFT they are also called cells.
Now
the nth harmonic can be expressed as n times f0.
this
is also equal to
Now
we rewrite the Fourier Transform substituting above expression for fn
in Eq. 2.
t’s cancel and we get,
-3-
The
above form of the Fourier Transform is called the Discrete Fourier Transform
(DFT). The Division by N is used to normalize the values.
DFT
is a special case of the Fourier Transform and is actually an approximation of
the real thing. The validity of the approximation is effected by the type of
waveform we are dealing with as well as the parameters fs and N.
Computing the DFT
The
process of computing the DFT is identical to computing the Fourier coefficients
we did in Tutorial 1.
Here
you need to know
1.
What is the sampling frequency of the target signal? Is the sampling frequency
large enough so that it covers all significant frequencies in the signal?
2.
How many samples do we need?
First
compute the fundamental frequency, and starting with the fundamental frequency
we multiply the discrete signal by a complex exponential and perform summation
on the result.
Do
you recall what it means to multiply by a complex exponential? How do you
interpret the following equation?
The
figure below shows what is happening in real-life.
Figure 6 - What we are really doing when we multiply by a complex
exponential
The
signal in fact is being split into two parts, 1. multiplied by a sine wave and
the other by a cosine of the same frequency. The resulting two signals are
orthogonal and are the result of multiplication with the complex exponential or
phasor.
DFT Step by step
Now
we compute the DFT of the signal in Fig 1.
Step
1 - Multiply the target signal in 7a by a cosine wave in 7b of frequency f0.
For this demonstration, we assume that f0 = 1. (Although only
cosines are shown, we do this for both sines and cosines and keep track of the
results separately.)
Figure 7a - The sampled target signal
Figure 7b - First Harmonic f1
= 1 * f0
Figure 7c - result of multiplying the first harmonic with the target
signal. The waveform has positive area.
The
multiplication gives us the waveform in 7cNow integrate this waveform over the
N samples. In a discrete case, we integrate by multiplying the sample amplitude
by the width of the base which is equal to t, the sample time, using the trapezoidal
rule. We are in effect adding up the areas of all the small gray rectangles in
Figure 7c. Each sample value is multiplied by t and these areas are summed.
The
result of the multiplication tells us something interesting. We see that the
resulting waveform is not even, so it has net area under it. This means that
there is a signal hiding in this frequency. What is the amplitude of this
frequency? That we know only when we complete the summation. The result of the summation
gives us the amplitude of this harmonic in the target signal.
Step
2: Now multiply the Signal in 7a with the second harmonic as shown in Figure
7d. The multiplication gives the waveform in Figure 7e.
Figure 7d - 2nd Harmonic
f2 = 2 * f0
Figure 7e - Result of multiplying the 2nd harmonic with the target
signal. The waveform has no area.
The
waveform of 7e is even, which means that the summation of the little gray
rectangles will give zero area. Since it has no net area means there is nothing
of interest here.
Let’s
go to the next harmonic. Now multiply the target signal with the 3rd
harmonic as in Fig 7f. The resulting waveform is shown in Fig 7g.
Figure 7f - 3nd Harmonic
f2 = 3 * f0
Figure 7g - Result of multiplying the 3rd harmonic with the target
signal. The waveform has positive area.
Once
again, we get a waveform that has a DC component, so we will do the summation
to figure out the amplitude.
In
Figure 7h and 7i we show the results of multiplication with harmonics 4 and 5.
For 5th harmonic we once again get a positive DC term. We will stop here because we already know
that our signal only has three frequencies. And we have found them all.
Figure 7h - Result of multiplying the 4th harmonic with the target signal.
The waveform has no area.
Figure 7i - Result of multiplying the 5th harmonic with the target
signal. The waveform has positive area.
Fast Fourier Transform (FFT)
The
DFT although clear and easy to compute requires a good many calculations. For
each harmonic, we have N+1 multiplications and we do this N times, giving us N2
+ N calculations. A 256 point DFT would require 65792 calculations. The
algorithm for Fourier Transform existed for more than 200 years before it came
into widespread use mostly because we could not cope with such large number of
computations. The algorithm waited for the development of the
microprocessor.
In
1948, Cooley and Tukey and came up with a computational breakthrough called the
Fast Fourier Transform algorithm. This was an ingenious manipulation of the
inherent symmetry of the calculations. It now allowed the computation of a N
point DFT as a function only 2N instead of the N2. So a 256 point
DFT would only require 512 calculations, a huge improvement from 65792 calculations
doing it the laborious way. The algorithm was quickly and widely adopted and is
the basis of all modern signal processing.
Most
DSP books spend a lot of time going through the mechanism of the FFT and how it
is computed. Huge butterfly figures in our books help to confuse and confound
us as to what is actually going on. Most of us just give up on it. Let me
alleviate your guilt. Although extremely important in itself, the understanding
of the mechanism of the FFT is not all that important. So, I am skipping over
the details of its implementation.
The
method has been programmed in all sort of software and we can safely skip it
without impacting our understanding of the application of the DFT and the FFT.
The
main thing one needs to know about the FFT is that it works only with sample
numbers that are powers of 2, such as 16, 32, 64 etc.. We cannot do a FFT on an
arbitrary number of samples as we can with DFT which is a generic process. The
FFT is a DFT with constraints on the number of samples.
The
other thing about the FFT process to know is that it allows zero-padding. Let’s
say we have 28 samples and we wish to do a DFT via the FFT, we can do two
things, 1. we can do a 16 point FFT and discard the remaining 12 points or we
can insert four zeros at the end so we have 32 points. Now we can do a 32 point
FFT. The zero-padding provides us better resolution but does not provide any
extra information. The frequency detected is still a function of the original N
samples and not the zero-padded length, although the FFT does look a lot better. And looks count.
FFT of a signal
Now
we will do an example using the signal
of Eq 1. This signal has three
frequency components. I have sampled it at the rate of 16 samples per second.
We will do several FFTs, from lengths 16 to 256 samples and examine the impact
of the sampling frequency and the number of samples.
Before
we do that let’s review the process of doing a FFT as shown in Fig. 8
1.
We
sample our continuos signal with an infinite length impulse train of frequency
fs.
2.
We
then limit the result to a certain number of samples by multiplying the result
of step 1 by a rectangular function of period T. All values outside this window
are now zero.
3.
We
realize that the FT of this limited signal is continuos because the Fourier
Transform of all periodic-discrete signals is continuos. We are after the
discrete Fourier Transform, so to do that we sample the limited signal in the
frequency domain.
4.
This
results in a replicated signal. We first windowed the signal to make it
non-periodic and then after windowing it, we have once again made it periodic.
This assumes that the windowed section is a true periodic representative of the
whole signal. And that by making it periodic again, we get the original signal
back. This is of course, just an assumption, which is rarely true for real
signals.
5.
The
Fourier Transform of this new periodic signal contains the effects of the
windowing function and as such is smeared. This smearing is most pronounced
when the window size is not an integral number of the signal period.
Figure 8 - DFT/FFT of a signal
Doing FFT of the signal g(t)
of Fig. 1
We
sample g(t) at the rate of 16 samples per second. Each frequency component of
g(t) is then also sampled at 16 samples per second. In each 16 samples, there
are an integer number of cycles of each of the three frequencies as shown in
Figure 9. (There is nothing magic about these numbers, all integers frequencies
will give an integer number of cycles per one cycle of the sampling frequency.
But had one of the frequencies been 3.5 instead of 3, we see than then we would
get a non integer number of cycles per 16 samples.)
fa = 1, contains 1 complete cycle in 16
samples
fb = 3, contains 3 complete cycles in 16
samples
fc = 5, contains 5 complete cycles in 16
samples
Figure 9 - Each
component of g(t) (fa), (fb) and ( fc) and the
composite signal g(t) which is the sum of (a), (b) and ( c) shown in (d) are
sampled at 16 times per second. When replicated the 16-sample pieces form the
original signal without discontinuities.
N
which is a multiple of 16 then always contains an integer number of cycles of
each of these frequencies.
In
the following examples, we show the effect of fs and N on the
results we obtain. Remember that we want the spectrum in Figure 2 based on
theory.
Case
1. We decide to use only 16 of the collected samples. Fs/N = 1 This
means that each frequency varies by 1 Hz. We see seven points. The values at
frequencies 1, 3 and 5 are correct but in between this is a terrible looking
representation of the real FT. All spectrums shown are one-sided and were done
in SPW.
Figure 10a - FFT of a three-frequency signal with N = 16
Case
2. Here w used 32 samples instead of 16. Now we have increased the resolution
to fs/N = 16/32 = .5. Now we get 14 points, each .5 Hz away.
Figure 10b - FFT of a three-frequency signal with N = 32
Case
3: N = 64 samples Now fs/N =
16/64 = .25. This Fourier Transform is
looking much better.
Figure 10c - FFT of a three-frequency signal with N = 64
Case
4: N = 128, The FFT looks nearly like the theory.
Figure 10d - FFT of a three-frequency signal with N = 128
Case
5: N = 256, This FFT looks quite satisfactory.
Figure 10e - FFT of a three-frequency signal with N = 256
What
conclusion can we draw from these? It is clear that the factor fs/N
has the largest impact. So we can always improve the FFT by increasing the size
of the FFT. The smaller the fs/N parameter, the better the
resolution of our spectrum.
However,
in this example all things were perfect. The window size contained exactly an
integer cycles of the target signal.
What
happens when the window size does not cover an integer multiple number of
cycles?
Case 2. Sampling frequency
not an integer multiple of the FFT length
In
this case, I have sampled the signal with fs = 20 instead of 16. We
will use 32 as the smallest FFT size. Let’s look at what portion of the signal
we capture in these 32 samples.
fa = 1, contains 1 complete cycles in 20
samples or 1.6 cycles per 32 samples
fb = 3, contains 3 complete cycles in 20
samples or 4.8 cycles per 32 samples
fc = 5, contains 5 complete cycles in 20
samples or 8 cycles per 32 samples
(please
note that the 3rd frequency does have an integral number of cycles
per 32 samples.)
Figure 11 - Each
component of g(t) is sampled 20 times per second. When replicated the
20-sample pieces do
not form the original signal. The replicated signal has discontinuities.
Now
we do FFT of this signal. N which is
always a multiple of 16 contains a non-integer number of cycles of the 1st
and the 2nd frequencies. We will see smearing in these two
frequencies caused by the lack of a truly periodic signal.
Figure 12a - FFT of a three-frequency signal with N = 32
Figure 12b - FFT of a three-frequency signal with N = 64
Figure 12c - FFT of a three-frequency signal with N = 128
Figure 12d - FFT of a three-frequency signal with N = 256
Figure 12e - FFT of a three-frequency signal with N = 512
Figure 12f - FFT of a three-frequency signal with N = 1024
Compare
this FFT to the 256 point FFT for the first case. Both signals are the same,
the only difference is the sampling frequency. In the second case the FFT
length no longer covers an integer number of cycles and so discontinuities are
present in the signal.
However,
it is interesting to note that the representation is not useless, and is in
fact quite good. The last frequency which we know is 5, gives us quite good
results. Why? Because frequency of 5 divides evenly into the FFT lengths of 32,
64, 128, 256, 512 and 1024 so it is represented faithfully.
A
general rule can be stated as:
where
fs = sampling frequency, fsig = signal frequency, 2n
are all numbers that are powers of 2 such as 16, 32, 64, 128, etc.
Example: fs = 20
fa
= 1
No
solution to above equation for N.
fb
= 3
No
solution to above equation for N.
fc
= 5
for
n = 5, we can solve the above equation. The minimum FFT size is 2n
or 32 and for this signal frequency and sampling frequency combination, we will
get a good FFT as we saw in the example FFT.
We
can also draw this conclusion: if we do a large enough FFT, we can to a large
extent overcome the problems associated with unknown truncations and can get
acceptable results. A concept called windowing the function with something
other than the rectangular function we used also helps to reduce these effects.
We will discuss this in the next issue.
Zero padding
In
most case, the FFT length is always less than the number of samples available
on hand. But sometimes we do an FFT that is larger that the number of samples.
In this case, the algorithm appends a bunch of zeros to the end of the signal
to artificially make it longer. Doing so sometimes gives better resolution for
the signal space and also gives us “good looking” spectrums. Whether or not to
zero pad is a subjective thing and depends on how much data you have and what
type of signal it is.
Figure 13a - FFT of a three-frequency signal with actual samples = 2000,
Zero-padded to N = 2048
Figure 13b - FFT of a three-frequency signal with actual samples = 2000,
Zero-padded to N = 4096
Figure 13c - FFT of a three-frequency signal with actual samples = 2000,
Zero-padded to N = 8192
This
is a good looking spectrum and demonstrates that the resolution, fs/N
is an important parameter of how good of a FFT we will get. The above FFT
created from 2000 points and then zero padded to 8192 is really very good. The
fact that the input signal was not a faithful reproduction of the actual signal
is nearly inconsequential.
How
wonderfully forgiving and idiot proof DFT/FFT is!
Next:
Periodogram, aliasing, leakage and windowing and other
approximations/modifications of the Fourier Transform
By:
Charan Langton
Complextoreal.com
mntcastle@earthlink.net